#include "Interpolation.h"

#include <cassert>
#include <stdint.h>

namespace NumericalAnalysis
{
    double CompositeTrapezoidRule(const double a, const double b, const double h, double (* const f) (const double))
    {
    	assert(b >= a && "CompositeTrapezoidRule: Wrong interval.");

    	uint32_t n = static_cast<uint32_t>( (b - a) / h );

    	double sum = 0.0;
    	for(uint32_t i = 1; i < n; ++i)
    		sum += f(a + static_cast<double>(i) * h);

    	return (h / 2.0) * (f(a) + 2.0 * sum + f(b));
    }

    double CompositeSimpsonRule(const double a, const double b, const double h, double (* const f) (const double))
    {
    	assert(b >= a && "CompositeTrapezoidRule: Wrong interval.");

    	uint32_t n = static_cast<uint32_t>( (b - a) / h );

    	double sum1 = 0.0;
    	for(uint32_t i = 2; i <= (n / 2); ++i)
    		sum1 += f(a + (2.0 * static_cast<double>(i) - 2.0) * h);

    	double sum2 = 0.0;
    	for(uint32_t i = 1; i <= (n / 2); ++i)
    		sum2 += f(a + (2.0 * static_cast<double>(i) - 1.0) * h );

    	return (h / 3.0) * (f(a) + 2.0 * sum1 + 4.0 * sum2 + f(b));
    }
}
